/*
**  File: AquaCommon.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht, the Netherlands and the BioMagResBank, 
**  University of Wisconsin-Madison, USA.

**
**  AUTHOR(S) : Jurgen F. Doreleijers
**  VERSION   : 3.2
**  DATE      : January 2001
**
**  This file contains miscellaneous functions for 
**     AquaAssign, AquaRedun, AquaSatur.
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"
#include "AquaTypes_Assign.h"

/****************  function prototypes  ****************/

#include "MenuFuncts.h"
#include "Functs.h"

#include "AquaCommon.h"
#include "nrutil.h"

/****************  globals  ****************/

extern Boolean  iaf;                   /* import from Functs */
extern int      terseness;             /* import from Functs */
extern Card     msg_string;            /* import from Functs */


/* BEGIN_EXTH */


/****************  Check_Redundant  ****************/
/*
**  Checks a list of DRestraint on redundancy.
*/


int Check_Redundant ( DRestraint *DRs, int count, int notkeep_count,
		Boolean detail,
		Boolean save_pseudos, Boolean del_red_l, Boolean del_red_u, 
		Boolean use_distance, float thresh )

{

#define RES_COUNT		    20
#define MAX_ATOMS_RES		    29
#define CHANGE_COUNT		    16
#define MIN_LOW                      0. 
#define MAX_UPP                     10.5 /* Highest found is 10.3 */
#define JFDBUG                       false

/* 
Compared  set of 5000  with  set of 10000 structures
Checked: 5590 out of 15264 distances
Average abs(difference) : 0.0013 +/- 0.0122
Average abs(difference) : 0.0270 +/- 0.2797 %
Maximum UP diff.:    0.2  and    4.0 %
Maximum LO diff.:   -0.1  and   -5.3 %
*/


    int       i, j, k, x;
    int       res_num, at_a, at_b, checked_count;
    
    Boolean   low_red, high_red, found, found1, found2,  skip;

    char      *PathName, Rdata[] = "redundant.data";
    Filnam    FileName;
    FILE      *fp;
    Card      txt[MAX_ATOMS_RES+1];             /* LAATSTE +1 voor ? */
    Boolean   PseudoA, PseudoB;
    Boolean   lowerbound, upperbound;
    double    percent, LowThres, UppThres;
    

    double Upp[RES_COUNT][MAX_ATOMS_RES][MAX_ATOMS_RES];
    double Low[RES_COUNT][MAX_ATOMS_RES][MAX_ATOMS_RES];

    ID    Residues[RES_COUNT];
    int Atom_Count[RES_COUNT];
    ID       Atoms[RES_COUNT][MAX_ATOMS_RES];
    
    ID ChangeAtoms[CHANGE_COUNT][3]= {
    
	{"ALA", "QB",  "MB"},   
	{"ASN", "QD2", "QD"},   
	{"GLN", "QE2", "QE"},   
	{"GLY", "HA1", "HA2"},   
	{"GLY", "HA2", "HA3"},   
	{"ILE", "QG1", "QG"},   
	{"ILE", "QG2", "MG"},   
	{"ILE", "QD1", "MD"},   
	{"LEU", "QD1", "MD1"},  
	{"LEU", "QD2", "MD2"},  
	{"LEU", "QQD", "QD"},   
	{"MET", "QE",  "ME"},   
	{"THR", "QG2", "MG"},   
	{"VAL", "QG1", "MG1"},  
	{"VAL", "QG2", "MG2"},  
	{"VAL", "QQG", "QG"}, };  
	
    int bad_count       = 0; 
    int inter_count     = 0; 
    int bad_resat_count = 0; 
    int pseudo_count    = 0; 
    int redundant_count = 0;
    int same_count      = 0;

    int both_red_count  = 0;
    int low_red_count   = 0;
    int high_red_count  = 0;

    int both_imp_count  = 0;
    int low_imp_count   = 0;
    int high_imp_count  = 0;

    int both_count      = 0;
    int low_count       = 0;
    int high_count      = 0;

    int notkeep_count_actual = 0;

    if ( terseness < 9 ) {
        printf("\nChecking the redundancy of intra-residual restraints\n\n");
    }

        
    PathName = GetVar( "AQUADATADIR", true );
    strcpy( FileName, PathName );
    strcat( FileName, DIRSEP );
    strcat( FileName, Rdata );
/* the total length of the string should be checked against size_FILNAM */
    F_OPEN_PROTECTED (fp,FileName,"r");

    if ( terseness < 6 ) {
        printf("Reading lower and upper limits from file: %s\n", FileName);
    }
		
    for (i=0;i<RES_COUNT;i++) {

	FlushLine( fp ); /* Blank line */
	
	k = 0;
	while ( ReadWord( txt[k], size_LINE, true, fp ) ) {
	    k++;
	}
	Blank( Residues[i], size_ID );
	CopyChars( Residues[i], txt[1], strlen(txt[1]) );

	sscanf   ( txt[3], "%d", &Atom_Count[i]);
	FlushLine( fp ); 
	
	k = 0;
	while ( ReadWord( txt[k], size_LINE, true, fp ) ) {
	    Blank( Atoms[i][k], size_ID );
	    CopyChars( Atoms[i][k], txt[k], strlen(txt[k]) );
	    k++;
	}
	if ( k != Atom_Count[i] ) {
	    printf("%d %d %d\n", i, k, Atom_Count[i] );
	    nrerror("ERROR : Bad read0");
	}
	FlushLine( fp ); 
	

	for (j=0;j<Atom_Count[i];j++) {

	    k = 0;
	    while ( ReadWord( txt[k], size_LINE, true, fp ) )
		k++;
	
	    FlushLine( fp );
	
	    if ( k != Atom_Count[i] ) {
		printf("%d %d %d\n", i, j, k);
		nrerror("ERROR : Bad read1");
	    }
	
	    for (k=0;k<Atom_Count[i];k++) {
		sscanf   ( txt[k], "%lf", &Low[i][j][k]);
		if ( ( Low[i][j][k] < MIN_LOW ) || ( Low[i][j][k] > MAX_UPP ) ) {
		    printf("%d %d %d\n", i, j, k);
		    nrerror("ERROR : Bad1 ");		    
		}
	    }
		
	    k = 0;
	    while ( ReadWord( txt[k], size_LINE, true, fp ) )
		k++;

	    FlushLine( fp );

	    if ( k != Atom_Count[i] ) {
		printf("%d %d %d\n", i, j, k);
		nrerror("ERROR : Bad read2");
	    }
    
	    for (k=0;k<Atom_Count[i];k++) {
		sscanf   ( txt[k], "%lf", &Upp[i][j][k]);
		if ( ( Upp[i][j][k] < MIN_LOW ) || ( Upp[i][j][k] > MAX_UPP ) )  {
		    printf("%d %d %d\n", i, j, k);
		    nrerror("ERROR : Bad2 ");
		}		    
		if (j < k )
		    if ( (Low[i][j][k] < TWO_PROTON_RADII) && 
		         (Upp[i][j][k] > TWO_PROTON_RADII) )
			Low[i][j][k] = TWO_PROTON_RADII;
	    }
	    	    
	    for (k=0;k<Atom_Count[i];k++)
		if (  Low[i][j][k] > Upp[i][j][k] )  {
		    printf("%8.3f %8.3f %d %d %d\n", Low[i][j][k], Upp[i][j][k], i, j, k);
		    nrerror("ERROR : Low > Upp ");		    
		}

	}
    }
    fclose(fp);

    /* Change MOLMOL NAMES */

    for (i=0;i<CHANGE_COUNT;i++) {
        found1 = false;
	j=0;
	while ( !found1 ) {
	    if ( memcmp (ChangeAtoms[i][0], Residues[j], strlen(ChangeAtoms[i][0]) ) == 0 ) {
		found1 = true;
		
		k=0;
		found2 = false;
		while ( !found2 ) {
		    if ( memcmp (ChangeAtoms[i][1], Atoms[j][k], strlen(ChangeAtoms[i][1])) == 0 ) {
			found2 = true;
			Blank( Atoms[j][k], size_ID);
			CopyChars( Atoms[j][k], ChangeAtoms[i][2], strlen(ChangeAtoms[i][1]) );
		    }
		    k++;
		    if ( k>Atom_Count[j] ) {
			PRINTID( ChangeAtoms[i][1] );
			nrerror("\nCouldn't find this atom");
		    }
		}
	    }
	    j++;
	    if (j>RES_COUNT) {
		PRINTID( ChangeAtoms[i][0] );
		nrerror("\nCouldn't find this residue");
	    }
	}
    }

    /* Symetrize */
    for (i=0;i<RES_COUNT;i++) {
	if ( JFDBUG ) {
	    printf("\nRESIDUE: ");
	    PRINTID( Residues[i] );  
	    printf("WITH %d ATOMS\n", Atom_Count[i]);
	    for (j=0;j<Atom_Count[i];j++)
		PRINTID( Atoms[i][j] );
	    printf("\n");
	}
	for (j=0;j<Atom_Count[i];j++) {
	    for (k=0;k<j;k++) {
		Upp[i][j][k] = Upp[i][k][j];
		Low[i][j][k] = Low[i][k][j];
	    }
	    if ( JFDBUG ) {
		for (k=0;k<Atom_Count[i];k++)
		    printf("%4.1f ", Low[i][j][k]);
		printf("\n");
		for (k=0;k<Atom_Count[i];k++)
		    printf("%4.1f ", Upp[i][j][k]);
		printf("\n");
	    }
	}
    }

    if ( JFDBUG ) 
	printf("\n\n");


   
    if ( terseness < 6 ) {
        printf("\n");
        printf("Bounds between atoms are checked with distances that are\n");
        printf("calculated from 10,000 more or less random structures\n"); 
        printf("using torsion angle dynamics in XPLOR.\n");
        printf("No (even not vdw) interactions were used in those calculations!\n");
        printf("The lower bound checking as a result is too generous,\n");
        printf("eg more lower bounds might be redundant because the lowest\n");
        printf("distances found is not always possible due to steric hindrance.\n");
        printf("Corrected this somewhat by using a minimum lower bound of:\n");
        printf("%3.1f Angstrom.\n\n", TWO_PROTON_RADII);

        printf("OPTIONS SET:\n");

        printf("    Threshold for redundancy = %.1f PERCENT.\n", thresh);
        printf("    detail is ");
        if ( detail ) {
	    printf("ON\n");
        } else {
	    printf("OFF\n");
        }
        printf("    Lower bounds <= %3.1f Angstrom are considered as not being a lowerbound.\n", 
	    TWO_PROTON_RADII);

        if ( del_red_u )
	    printf("    Will consider restraints redundant, even if lower bound is not redundant.\n");	

        if ( del_red_l )
	    printf("    Will consider restraints redundant, even if upper bound is not redundant.\n");	

        if ( save_pseudos )
	    printf("    Will not delete non-redundant restraints with one or two pseudo atom(s).\n");

        if ( use_distance )
	    printf("    Will use the distance value if no (real) lowerbound is given.\n");	
    }

    if ( detail && (terseness < 9) ) {
	printf("\n\nCh. Resid.  Atom1 Atom2  Possiblebounds Restraintbound(s) Correction\n");
    }
    
    
    for (x=0;x<count;x++) {

        res_num = at_a = at_b = -1;

	if ( DRs[x].keep == false ) {
	    notkeep_count_actual++;
	    continue;
	}

	lowerbound = ( (DRs[x].stat == bnd_lower) || (DRs[x].stat == bnd_uplo) 
	                                          || (DRs[x].stat == bnd_uld ) );
	upperbound = ( (DRs[x].stat == bnd_upper) || (DRs[x].stat == bnd_uplo) 
	                                          || (DRs[x].stat == bnd_uld ) );

	if ( lowerbound ) {
	    if ( DRs[x].lo_bnd <= TWO_PROTON_RADII )
		lowerbound = false;	
	    if ( !lowerbound && (DRs[x].stat == bnd_uld ) && use_distance ) {
		lowerbound = true;
		DRs[x].lo_bnd = DRs[x].dist;
	    }
	}
    
	if ( lowerbound && upperbound && (DRs[x].lo_bnd > DRs[x].up_bnd) ) {
	    printf("ERROR : lo_bnd > up_bnd");   
	    printf(" BAD Restraint # %3d :low dist up= %6.3f %6.3f %6.3f\n", x, 
	    DRs[x].lo_bnd, DRs[x].dist, DRs[x].up_bnd);
	    bad_count++;
	    continue;
	}

    	if ( ( memcmp ( DRs[x].atom_A.chain, DRs[x].atom_B.chain, size_ID ) != 0 ) ||
	     ( DRs[x].atom_A.res_num_orig !=  DRs[x].atom_B.res_num_orig ) ) {
		inter_count++;
		continue;
	}

	
	if ( detail && (terseness < 9) ) {	
	    PRINTIDv( DRs[x].atom_A.chain    ); printf(" ");		 	
	    printf("%3d ", DRs[x].atom_A.res_num_orig);
	    PRINTID( DRs[x].atom_A.res_nam  ); printf(" ");		 	
	    PRINTID( DRs[x].atom_A.atom_int ); printf(" ");	 	
	    PRINTID( DRs[x].atom_B.atom_int ); printf(" ");	 	
	}

        /* Check if residue and atom names can be found */
	
	found = false;
	skip  = false;
	
	for (j=0;j<RES_COUNT;j++) {
	    if ( memcmp ( Residues[j], DRs[x].atom_A.res_nam, size_ID ) == 0 ) {
		found = true;
		res_num = j;
	    }
	}
	if ( ! found )  {
	    skip = true;
	    if ( detail && (terseness < 9)  ) {	
		printf("Couldn't find Residue ");
	    }
	}
	
	if ( !skip ) {
	    found = false;
	    for (j=0;j<Atom_Count[res_num];j++) {
		if ( memcmp ( Atoms[res_num][j], DRs[x].atom_A.atom_int, size_ID ) == 0 ) {
		    found = true;
		    at_a = j;
		}
	    }
	    if ( ! found )  {
		skip = true;
		if ( detail && (terseness < 9)  ) {	
		    printf("Couldn't find Atom_A ");
		}
	    }
	    found = false;
	    for (j=0;j<Atom_Count[res_num];j++) {
		if ( memcmp ( Atoms[res_num][j], DRs[x].atom_B.atom_int, size_ID ) == 0 ) {
		    found = true;
		    at_b = j;
		}
	    }
	    if ( ! found ) {
		skip = true;
		if ( detail && (terseness < 9)  ) {	
		    printf("Couldn't find Atom_B ");
		}
	    }
	}
	
	if ( skip ) {
	    if ( detail && (terseness < 9) ) {	
		printf("BAD RESIDUE OR ATOM\n");
	    }
	    bad_resat_count++;
	    continue;
	}	

	high_red = true;
	low_red  = true;

	if ( detail && (terseness < 9) ) {
	    printf("lo %4.1f ", Low[res_num][at_a][at_b]);
	    printf("hi %4.1f ", Upp[res_num][at_a][at_b]);
	}
					
	LowThres = Low[res_num][at_a][at_b]*thresh/100.0;
	UppThres = Upp[res_num][at_a][at_b]*thresh/100.0;

	if ( lowerbound ) {
	    if ( detail && (terseness < 9) ) {
		printf("LO %5.2f ", DRs[x].lo_bnd);
	    }
	    if ( DRs[x].lo_bnd > Low[res_num][at_a][at_b] - LowThres )
		low_red  = false ;
	}
	    
	if ( upperbound ) {
	    if ( detail && (terseness < 9) ) {
		printf("HI %5.2f ", DRs[x].up_bnd);
	    }
	    if ( DRs[x].up_bnd < Upp[res_num][at_a][at_b] + UppThres )
		high_red = false;
	} 
	
	    	
		
	if ( save_pseudos ) {
	    PseudoA = false;
	    PseudoB = false;

	    if ( memchr( Atoms[res_num][at_a], 'M', size_ID) != NULL )
		PseudoA = true;
	    if ( memchr( Atoms[res_num][at_a], 'Q', size_ID) != NULL )
		PseudoA = true;
	    if ( memchr( Atoms[res_num][at_b], 'M', size_ID) != NULL )
		PseudoB = true;
	    if ( memchr( Atoms[res_num][at_b], 'Q', size_ID) != NULL )
		PseudoB = true;
		
	    if ( PseudoA || PseudoB  ) {
		if ( detail && (terseness < 9) ) {
		    printf(" PSEUDO SKIPPED\n"); 
		}
		pseudo_count++;
		continue;
	    }
	}


	if ( lowerbound ) {
	    low_count++;
	    if ( upperbound ) 
		    both_count++;
	}
	if ( upperbound ) 
	    high_count++;


	if ( Low[res_num][at_a][at_b] ==  Upp[res_num][at_a][at_b] ) {
	    if ( detail && (terseness < 9)  ) {
		printf("NO FREEDOM THUS BOTH BOUNDS REDUNDANT\n");
	    }
	    DRs[x].keep = false;
	    same_count++;
	    redundant_count++;
	    if ( lowerbound ) {
		low_red_count++;
		if ( upperbound ) 
			both_red_count++;
	    }
	    if ( upperbound ) 
		high_red_count++;
	    continue;
	}
		
	    
	if ( lowerbound && low_red ) {
	    if ( upperbound && high_red ) {
		if ( detail && (terseness < 9)  ) { 
		    printf("BOTH REDUNDANT ");
		}
		high_red_count++;
		low_red_count++;
		both_red_count++;
	    } else {
		if ( detail && (terseness < 9)  ) { 
		    printf("LOWER REDUNDANT ");
		}
		low_red_count++;		
	    }
	}
	if (( upperbound && high_red ) && !(lowerbound && low_red) ) {
	    if ( detail && (terseness < 9)  ) {
		printf("UPPER REDUNDANT ");
	    }
	    high_red_count++;
	}
	
	if ( lowerbound && ( DRs[x].lo_bnd > Upp[res_num][at_a][at_b] + UppThres ) ) {
	    if ( upperbound && ( DRs[x].up_bnd < Low[res_num][at_a][at_b] - LowThres ) ) {
		if ( detail && (terseness < 9)  ) { 
		    printf("BOTH ARE IMPOSSIBLE ");
		}
		high_imp_count++;
		low_imp_count++;
		both_imp_count++;
	    } else {
		if ( detail && (terseness < 9)  ) { 
		    printf("LOW=IMPOSSIBLE ");
		}
		low_imp_count++;
	    }
	}
	if ( upperbound && ( DRs[x].up_bnd < Low[res_num][at_a][at_b] - LowThres ) 
	     && ! ( lowerbound && ( DRs[x].lo_bnd > Upp[res_num][at_a][at_b] + UppThres ) ) ) {
		if ( detail && (terseness < 9)  ) {
		    printf("HIGH=IMPOSSIBLE ");
		}
		high_imp_count++;
	}	
	
	if ( (  low_red  &&  high_red  ) || ( high_red && del_red_u ) || ( low_red && del_red_l ) ) {
	    redundant_count++;
	    DRs[x].keep = false;		
	}
	if ( detail && (terseness < 9)  ) {
	    printf("\n");
	}
    }
    
    
    if ( notkeep_count != notkeep_count_actual ) {
	printf("Found %d and %d restraints not to keep, doesn't match\n",
	    notkeep_count, notkeep_count_actual );
	nrerror("\nError in Check_Redundant");
    }
    
    printf("Total of %d restraints scanned.\n", count );

    checked_count = count - inter_count - pseudo_count 
                    - bad_resat_count - bad_count - notkeep_count_actual;

    if ( (double) (checked_count) == 0 )
	percent = 0;
    else
	percent = 100.0 *  redundant_count / (double) (checked_count);

    printf("CHECKED  %4d INTRA residual restraints.\n", checked_count);
    printf("Found         %4d REDUNDANT restraints, %7.2f percent\n", 
        redundant_count, percent);

    printf("                    %4d with lower bound\n", low_count);
    printf("                    %4d with upper bound\n", high_count);
    printf("                    %4d with both bounds\n", both_count);
    printf("                    ****\n");

    printf("                    %4d with redundant lower bound\n", low_red_count);
    printf("                    %4d with redundant upper bound\n", high_red_count);
    printf("                    %4d with both bounds redundant\n", both_red_count);
    printf("                    ****\n");

    printf("                    %4d with impossible lower bound\n", low_imp_count);
    printf("                    %4d with impossible upper bound\n", high_imp_count);
    printf("                    %4d with both bounds impossible\n", both_imp_count);
    printf("                    ****\n");

    printf("                    %4d with same read Upp & Low\n", same_count);

    printf("Skipped  %4d \"not to keep\"  restraints.\n", notkeep_count_actual);
    printf("Skipped  %4d INTER residual restraints.\n", inter_count);
    printf("Skipped  %4d INTRA residual restraints with pseudo atom(s).\n", pseudo_count);
    printf("Skipped  %4d unknown INTRA residual restraints.\n", bad_resat_count);
    printf("Skipped  %4d restraints with lower bound > upper bound.\n", bad_count);

    return ( redundant_count );
}



void Calculate_Corrections (Extra_Distinfl *extra_distinfl, 
				Distinfl *dist_infl, int dist_count)
{
/*  Based on the type of pseudo atom; as defined in distinfl, from stereoLIB;
    set a correction in the extra_distinfl list for both pairs individually.

$ REMARK: 	   Only pseudo's of the type 1,3,4, and 5 can be used.
$ REMARK: 		0 undefined  (*)
$ REMARK: 		1 CH2 or NH2 (2) 		2 methyl or NH3 (3)
$ REMARK: 		3 double CH2/NH2 (4)		4 double methyl (6)
$ REMARK: 		5 Aromat with 2 H (2) 		6 Aromat with 4 H (4)

    Values taken from distance of pseudo atom to constituent atom as
    measured in INSIGHTII with ideal AA (leucine, tyrosine and arginine)
    from library.
    Values are rounded to higher value (eg 0.83->0.9) to be generous.
    The ones not to use are given a very high value (defined: BAD_CORRECT).
    A check is made to make sure these do not occur in both members.
    Type 6 was added, I have no idea why I left it out originally.
    The StereoLIB-aqua is also updated in version 3 with QR.
*/
    
    int i;
    double cor[7] = {BAD_CORRECT, 0.9, BAD_CORRECT, 1.5, 1.6, 2.2, 1.3 };
    
    for (i=0;i<dist_count;i++) {
    
	if ( ( Stereo_Assignable_Atom_Type( dist_infl[i].atype_A ) ) &&
	     ( dist_infl[i].count_A < 3 ) )
	    nrerror("ERROR in routine Calculate_Corrections; wrong atype A" );
			    
	if ( ( Stereo_Assignable_Atom_Type( dist_infl[i].atype_B ) ) &&
	     ( dist_infl[i].count_B < 3 ) )
		 nrerror("ERROR in routine Calculate_Corrections; wrong atype B" );

	if ( dist_infl[i].count_A == 1 ) 
	    extra_distinfl[i].cor_A = BAD_CORRECT; /* Should never be used */
	else 
	    extra_distinfl[i].cor_A = cor[ dist_infl[i].atype_A ];
				    
	if ( dist_infl[i].count_B == 1 ) 
	    extra_distinfl[i].cor_B = BAD_CORRECT; /* Should never be used */
	else 
	    extra_distinfl[i].cor_B = cor[ dist_infl[i].atype_B ];
	
    }
}		 



Boolean Stereo_Assignable_Atom_Type(int type)
{

    if ( (type == 0) || (type == 2) || (type == 6) )
	return (false);
    else {
	if ( (type < 0) || (type > 6) )
	    nrerror("ERROR in function Stereo_Assignable_Atom_Type" );
	return (true);
    }

}

